BI-TP 2011/44 



Thermal right-handed neutrino production rate 
in the non-relativistic regime 

M. Laine, Y. Schroder 

0J Faculty of Physics, University of Bielefeld, D-33501 Bielefeld, Germany 

D 

Abstract 

We consider the next-to-leading order thermal production rate of heavy right-handed neu- 
trinos in the non-relativistic regime mtop ^ ttT <^ M, where mtop refers to the electroweak 
scale. Rephrasing the problem in an OPE language and making use of different techniques 
than a previous analysis by Salvio et al, we confirm the general structure of their result and 
many of the coefficients. We also extend the analysis to the next order in the non-relativistic 
expansion, thereby revealing the leading non-trivial momentum dependence, as well as to 
NNLO in couplings, revealing the leading sensitivity to thermal resummations. Our results 
are expressed as a sum of simple "master" structures, which renders them a suitable starting 
point for determining the next-to-leading order rate also in the relativistic regime ttT ~ M. 
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1. Introduction 



Weakly interacting particles produced by scatterings taking place in a hot plasma could 
conceivably play a role in cosmology. For instance, some Dark Matter candidates, such as 
axions, axinos, gravitinos, or sterile neutrinos, could originate this way (for reviews see e.g. 
refs. [Il[2]). Another example is the so-called Leptogenesis model for explaining the Baryon 
Asymmetry of the Universe ^ , in which thermally produced out-of-equilibrium right-handed 
neutrinos could act as one of the main building blocks. 

Although the Leptogenesis scenario has been thoroughly studied (for reviews, see e.g. 
refs. [HE]), there are some ingredients in the analysis that do not appear to be on a sound 
theoretical footing. In particular, it has recently been pointed out that in the "ultrarelativis- 
tic" regime, meaning temperatures much higher than the mass of the right-handed neutrinos, 
vrT ^ M >mtop, the expressions that had been used are not correct even at leading order 
in the coupling constants (because of infrared sensitivity infinitely many loop orders con- 
tribute to the leading-order result in this kinematic regime), and that they had therefore 
underestimated the production rate [6]. On the other hand, in another recent contribution, 
it was shown that in the "non-relativistic regime" , mtop ^ ttT <^ M, the significance of loop 
effects had probably been overestimated in the literature, because they had been only partly 
accounted for, whereby an important cancellation was missed [7]. (More complicated issues 
like CP-violation could contain even larger uncertainties [8]-|13j.) 

In view of the mentioned developments, it might be useful to revisit the analysis of the right- 
handed neutrino production rate also in the "relativistic" regime, mtop ^ t^T ~ M, checking 
the importance of loop corrections through an explicit computation. Here, we take a first 
step in this direction, by formulating the setup of a next-to-leading order (NLO) analysis, 
and by carrying it out in the non-relativistic regime previously considered in ref. [7]. Given 
that the analysis of ref. [7] is exceedingly complicated, our goal was to reproduce its results by 
simpler methods, and also to extend them up to a higher order in an expansion in (vrT/M)^. 

It is important to understand that when discussing a non-relativistic regime, we place our- 
selves in the position of the particles produced, with an invariant mass K? = ^ (vrT)^, 
a kinetic energy k'^/{2M) ~ vrT, and therefore an average velocity k/M ~ ^JiiTlM <^ 1. 
From the plasma perspective, however, the four-momentum fC is an external probe, which 
has no dynamical effect. It simply represents a "hard" external scale, with , k,\k^ ± k\ ^ 
ttT. a proper tool for addressing this kinematic situation is the Operator Product Ex- 
pansion (OPE) [H], formulated within the context of thermal field theory only rather re- 
cently p2^J, and subsequently applied e.g. to correlators of the energy-momentum tensor of a 
Yang- Mills plasma [16\ I17j. Understanding the computation in this language has significant 
conceptual benefits, for instance by leading to a general framework for analyzing the infrared 
(IR) sensitivity of the results [K], a topic that remains otherwise to be studied empirically 
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2. Setup 



We consider a plasma made of Standard Model particles, in thermal equilibrium at a tempera- 
ture T, and interacting with right-handed neutrinos, of mass M, through Yukawa interactions, 
parametrized by a coupling hy . (It would be natural to assume the Yukawa couplings to build 
a 3 X 3 matrix, but our results are flavour-blind so we can simplify the notation without a loss 
of generality.) As long as the density of the right-handed neutrinos is below the equilibrium 
value, their production rate can be computed from first principles using a linear response 
or Kubo type analysis. Then, essentially, the production rate is determined by the cut (or 
imaginary part) of the self-energy of the right-handed neutrinos. 

An explicit derivation of the production rate has been presented in ref. Letting 
4) = icr2(j)*; denoting by £ a lepton doublet; and by a^,a^ the left and right projectors, 
Ol = (1 — 75)/2, Ojj = (1 + 75)/2, the differential production rate reads 



d7V(/C) \h 



UB I 



|2 



d^Xd^k (2^)3fcoTr{^a,[np(A;°-;u,)p(/C)+np(A;° + ^^)p(-/C)]aH,} + 0(|/i,B|') , (2.1) 

where h^^ is the bare neutrino Yukawa coupling; np{k^) = l/[exp(A;'^/T) + 1] is the Fermi 
distribution; /i^ denotes a leptonic chemical potential; p is the spectral function related to 
the composite operator (f)H; X = {x^,x.); and /C = {k^,k) is an on-shell four-momentum, 
with = + Af2. If the plasma is charge-symmetric, i.e. /i^ = 0, then p is symmetric in 
/C — )• — /C and the two terms can be combined; we make this assumption in the following. 
For a practical computation, it is convenient to employ Euclidean conventions. We define 



T. 



(2.2) 



where now X = (t, x); the vector K = (/c„,k) is Euclidean, with kn = (2n + 1)ttT, n G Z; 
and {■■■)rp refers to a thermal expectation value. Defining a differential "decay rate", F, in 
accordance with ref. [7] as 

diV(/C) _ 2np(fc0) 

dWk = 72^ ' ^^-^^ 

then eqs. (12. ip . (12. 2p together with the standard relation between the Euclidean correlator 
and the spectral function (p = ImllR = lm{IlE} k„^-i[kO+io+], cf. e.g. ref. [20]) imply that 

r(^) = ^im{n^(i^)},„^_,[,o+,o+]- (2-4) 

Here Im refers to a discontinuity, or cut, across the real A:'^-axis, and the relation is valid even 
on the non-perturbative level. 



^The notation in ref. [18] assumed a broken electroweak symmetry, but the derivation goes through also 
with a dynamical doublet. The derivation was extended to /if 7^ in ref. [19) . 
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To compute the correlator of eq. (j2.2p . let us establish conventions for the relevant part of 
the Standard Model. The Higgs field interacts according to the Euclidean Lagrangian 



+ htB q a^cj) t + h^j^t (p'^a^q 



(2.5) 



The mass parameter, m^, will mostly be omitted because we assume that vrT>mtop and, 
as already alluded to above, the NLO results turn out to be IR safe in this regime (cf. 
sec. [8]). Technically the calculation is performed as if we were in the symmetry restored 
phase although, as will become clear later on, this assumption can be modestly relaxed. 
When acting on the Higgs the covariant derivative takes the form 

[^M = {Smnd^ " ^52bT^„^^ + ^sTLi?^) , (2.6) 

where m,n £ {1,2}; are the SUl(2) and Uy(1) gauge fields, respectively; T° are 

Hermitean generators of SUl(2), normalized as Tr [T'^T''] = ^5"''; = ^Smn', and the 
couplings appearing are bare ones. The quark Yukawa interaction in eq. (|2.5p contains the 
doublet q = {t b)'^ . The leptons interact according to the Euclidean Lagrangian 



(2.7) 



where £ = {v ey^ . When acting on the lepton doublet, the covariant derivative reads 



(21 



3. Issues with regularization 

A loop computation in quantum field theory necessitates regularization, and by far the most 
convenient choice for this is the dimensional one, which we also adopt here. Unfortunately, 
with chiral gauge theories dimensional regularization leads to inevitable problems. Although 
well-known, we briefly remark on some of the issues in this section. 

In D dimensions, the normal (Euclidean) Dirac matrices can be taken to satisfy 

{7^,74 = 2g^, EE 2V , 7l = li., Tr(l)=4. (3.1) 
According to the 't Hooft - Veltman [21\ [22] convention a Hermitean 75 can be defined as 

75 = 70717273 , 75 = 1 • (3-2) 
An important implication from here is that |23[ 

aH,7/.aL = , (3.3) 
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where 

• (3.4) 

In the (free) kinetic part of eq. (12. 7p we might then imagine defining a Lagrangian as ij^dfj^i 
i^^d^i = ia^^^j^d^a^l + la^j^d^a^i, respecting gauge symmetry; but then the propagator is 
four-dimensional and not properly regularized. Or, we could keep the regularized form, 

= ^a^liMd/^aJ + la^-ffidf^a^i + ia^-ff,d^aj + la^"i^,d^a^i , (3.5) 

leading to the usual L'-dimensional propagator; but then gauge symmetry is (slightly) broken 
by the coupling of the different chiralities. Though extremely naive, these remarks already 
illustrate some of the complications encountered. Still, consistent computations are possible 
if counterterms and operator mixing are properly accounted for |23tl24j. and a closely related 
practical recipe has also been put forward [25] . 

In the prescription of ref . [25] , one can conveniently work with D-dimensional Dirac matri- 
ces after defining axial- vector currents with the structure 7/^75 — )• \ [7/^,75] = gr^^ivpo- 7i/7p7tT, 
pulling the antisymmetric tensors e outside the actual integrals and utilizing only their anti- 
symmetry inside the traces. The latter is made explicit by writing 

^[7^,75] = -^£tivpa{lvlpla-lalplv) ■ (3.6) 



As a further simplification, it has been argued in ref. [25] that a naively anticommuting 75 in 
combination with 7I = 1 can be used in fermionic traces containing more than one 75, with 
the exception of closed fermion loops. The consistency of this prescription has been verified 
up to 3-loop order in connection with singlet as well as non-singlet axial-vector operators in 
QCD [25^ 126]. We will employ this recipe as a "minimal crosscheck", commenting also on 
how it differs from the 't Hooft - Veltman scheme when closed fermion loops are present. 

It is important to realize, however, that in the computation of the spectral function accord- 
ing to eq. (|2.4p some of the ambiguities affecting the Euclidean correlator drop out. Indeed, 
a non-zero cut arises from logarithms or, in the context of dimensional regularization, from 
terms containing 1/e-poles (cf. eq. ()B.3P ). So, in terms of the Euclidean correlator, it is 
only necessary to get the 1/e-poles correct. Of course, the bare parameters also bring along 
1/e-poles, so lower-order graphs need to be worked out to a higher depth in the e-expansion. 



4. Leading order at zero temperature 

In order to get going, we start by an almost trivial step, performing a leading order analysis 
at zero temperature. This has the lucky feature of showing that at leading order, we are free 
from the subtleties of sec. [3] to all orders in e. 
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The starting point is to carry out Wick contractions in eq. (|2.2p . Then, although eq. (j3.3|) 
impUes that the external momentum appears as 



aR-^«L = -^aL, (4.1) 

rotational symmetry guarantees that we are free to choose K to have at most min(D — 1,3) 
non-zero spatial components. With this choice, Iji a^= Jji a^-, and 

Tr a^{i$ - if = 2k ■ {P - K) = 2K ■ {P - K) . (4.2) 

So we have the same expression as in Naive Dimensional Regularization (NDR), in which it 
is assumed, despite algebraic inconsistencies, that 75 anticommutes with all 7^. 



With eq. (j4.2p at hand, a few steps lead to 

^ P^{P-Kf ~ ^I'^'^^l'^^ [(P-K)2 " P2 - p2(p_ ^)2_ 



(4.3) 



where the additional factor 2 comes from the isospin trace, and we completed squares in the 
numerator. At zero temperature, with — )• Jp, only the last term contributes, and we 
obtair§ 

<)(i^) = -2|/..3pi^^^0 + ln|^ + 2 + O(e)) . (4.4) 
Subsequently, eq. (|2.4p yields 

ri^\lC) = \!^ + Oie), (4.5) 

where we have denoted /C^ = (k^)'^ — k'^; and the Euclidean four-momentum transforms as 
K'^ = k'^ + k'^ ^ -K? - i sign(/c°) 0+. The result in eq. (|4.5p agrees with eq. (3) of ref. |7j. 



For later reference, let us verify eq. (14. 3p using the prescription introduced below eq. (13. 6p . 
Performing the isospin trace, but keeping the Dirac trace for the moment, we get 

UEiK) = -2|/^.BpTr[7^fl,7.aJ^ ^pJi^Zxy ' ^^"^^ 

Now, noting that Tr (7^.aL7^aj^) = 5Tr{7^(7;, -757^^75) + l^l[lu,l5\}, the commutator part 
leads to ^ Eyp^K Tr [7^7p7cr7K— 7/x7k7o-7p] and vanishes once the trace is expressed as a product 
of metric tensors. In the first part, the rule of ref. [25] with two 75's leads immediately to 
[TM^LTf'^R] ~^ '^9fj.y, agreeing with eq. (|4.3p . In contrast, a strict use of eq. (13. 2p with two 
75's leads to the middle equation in eq. (|4.2p and necessitates a further argument. 



^The MS scale parameter is introduced in a usual way, inserting 1 = /i 2e^2<; cxpi^j-Ei) ^ 
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5. General method at non-zero temperature 



Proceeding to NLO, we illustrate some details by computing explicitly the "Higgs correction" , 
i.e. terms proportional to A. 



5.1. Wick contractions 

Carrying out the Wick contractions and making use of eq. ()4.2I) we obtain, in analogy with 
eq. gD, 



n 



(A) 



12|^,,n| Ap 



K-{K-P) 
K^ + {K- P f - P2 



P,Q 



Q^p^{K - py 



12\K 



1 



1 



+ 



p4 p2{^K-PY P^{K-PY 



(5.1) 



The factorized structures have well-known expressions in dimensional regularization |27j : 



Q 



1 

1 

P4 



-2e 



q H 

1 - 2e 



12 



1 + 2e 



In 



47rr 



-2e 



- + 21n(^ 
e V 47rr 



+ 1 + 



7e 



C(-i) 

+ 0(e 



+ 0(6')^ (5.2) 
(5.3) 



where n^{p) = l/[exp{p/T) — 1] denotes the Bose distribution, and a partial integration was 
carried out for obtaining the first representation of the latter term. (It is worth noting that 
this term has a power-like IR divergence which is not visible in dimensional regularization; 
the 1/e has an ultraviolet (UV) origin. Nevertheless, as long as dimensional regularization is 
applied consistently, there is no reason to worry; we return to this in the paragraph following 
eq. (|5.15p . and a comprehensive analysis is presented in sec. [H) However, the product of 
eqs. (j5.2p . (jS.Sp is i^-independent, so there is no cut, and no contribution to r(/C). 

5.2. Single-pole cut 

Non-trivial cuts arise from the 2nd and 3rd terms of eq. (j5.ip . Starting with the second term, 
the Matsubara sum can be carried out exactly, yielding 



1 



P^{K-PY JpP'^{K-Pf 



1 



+ 



V 



1 



{K - P)2 



(5.4) 



where the notation [...]p = \ Z]p„=±jp(---) corresponds to setting the thermal line on-shell; 
and the integration variable was renamed in the fermionic cut. The vacuum term has the 
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familiar structure of eq. (j4.4p . with a cut j^, whereas the thermal part can formally be 
expanded in p, given that the p-integration is exponentially convergent due to the thermal 
distributions: 



1 




1 




XK-py 


p 


K"^ -2K-P 


p 



1 2K-P A{K-Pf 



(5.5) 



This yields an OPE-type expansion. However, as long as we stay away from the light-cone 
-i[kO+iO+] 0)) this expansion has no cut. Therefore, a contribution only arises 
from the vacuum part of eq. (|5.4p , multiplying eq. (j5.2p : 



Im 



P,Q 



Q^P'^iK - Py ) k„^-i[ko+iQ+] 



T67rT2 



+ 0{e) 



(5.6) 



It is useful for later reference to repeat the analysis of the thermal part in a more detailed 
way. For this purpose we set k = and work out the cut exactly. Given that for A: = 



1 



1 



P;k=0 



{K-py 

we get 

'K[(A--P)^ 

So, for = Vfc2 + M2 > 0, 



1 



+ 



1 



1 



1 



1 



2kn \ kn — 2ip kn + 2ip J Aip \ kn — 2ip kn + 2ip 



(5.7) 



1 



P;fc=OJ k„^-i[kO+iO+] 



4p 



6{k^ - 2p) - S{k° + 2p) 



(5.8) 



Im 



UBjp) - nF{p) 
p 



1 



{K - Pf 



P:k=0 



1 f°° 

— J dp [nji{p) - np{p)]6{k^ - 2p) 



1 

16^ 



Tip 



A;0 



(5.9) 



In other words we do find a correction to the appearing in eq. (j5.6p but it is exponentially 
small, suppressed by e"'^"/^ in the regime A;° » T that we are interested in. 



5.3. Double-pole cut 

It remains to work out the 3rd term of eq. ()5.ip . This is slightly less trivial than the 2nd 
term because of the appearance of the double pole. Nevertheless the basic point is the same: 
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denoting Ep = p, E^p = |k — p|, we get 
1 



d 

dEl 



_d 
dEl) 



1 



{pl + El)[{kn-PnY + El^] 
1 



27r {pl + El)[{kn-po? + E^ 



+ 



nB{Ep) 

En 



{K-PnY + El 



L 



kp 



n-ejEkp) 



1 



pI + e^ 



K-P 



P P^{K-PY 



+ 



d n^iEp) ^ np{Ep) d 



dEl Ep 



{kr,-PnY + El^ 



(5.10) 



where we substituted p — )• k — p in the last term. The thermal part is more comphcated 
than before, but it is clear that it can again be expanded in inverse powers of l/K"^, and all 
resulting integrals are exponentially convergent in the UV (in the IR there are divergences 
like in eq. (j5.3p but this only leads to logarithms of T, not of K). Like in eq. ()5.5p . these 
terms give no cuts. The zero-temperature term, on the other hand, evaluates in dimensional 
regularization to 



p p^{K -py 



K-^'T{l + e)T^{l-e) 
eK^ (47r)2-^r(l - 2e) 



-2e 



Here 1/e is an IR divergence, but the cut is IR finite, and can be extracted like in eq. (14.4 
In total, then. 



Im 



P,Q 



Q2p4(^_p)2 



-i[kO+iO+] 



1 



+ 0(e) . 



(5.12) 



For the benefit of a reader nerved by the free use of dimensional regularization for handling 
IR divergences, let us repeat the derivation of eq. (|5.12|) in a more pedantic way, directly in 

= 4, avoiding IR divergences. Setting again k = for simplicity, and denoting Eq = 
p^ + ttiq, where tuq > is an IR regulator of the scalar line, we can rewrite eq. (jS.lOp as 



P^{K -Pf 



d 



lim „ , 



d 

lim - — 2 



1 



1 



1 



+ 



ikn+p + Eq 
1 



(pl+Ei)[{kn-pnr+p'] 
1 1 r 

1 - n^ip) + nji{Eo) 
n^{p) + n^iEo)]] . (5.13) 



+ 



1 

-ikn +P + Eo 
1 



Setting kn 



ikn+p-Eo -ikn+p-Eo_ 
-i[k^ + iO'^] and taking the cut, the denominators turn into 5-constraints like 
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in eq. (jS.Sp . If we choose mo < M, only one of them can get reahzed: 

.2 



Im 



P^{K - P) 

= lim / lh_s^-k° +p + Eo)\l-7iAp) + n^{Eo] 

mo->-o dniQ Jp 4:pEo I 



-i[fcO+iO+] 
.2 



lim 

mo^o drriQ I 167r 



2A;o ^ V 

Takmg the derivative we do reproduce eq. (|5.12|) . up to exponentiahy small terms. 

5.4. Putting everything together 

Inserting finally eqs. (15. 6p . (j5.12p into eq. (jS.lh . we obtain 



(5.14) 



>-i[fc"+iO+] \ Svr / 

This corresponds to the Higgs contribution in eqs. (16. 9|) . (l7.2p PI 

We now return to the issue of IR divergences. As has been stressed in the introduction, from 
the point of view of the thermal medium the right-handed neutrino acts as a "hard probe" ; 
any thermal effects on its production rate can be understood in the OPE language. But in the 
OPE language thermal effects amount to gauge-invariant condensates developing expectation 
values (15| . These, in turn, obtain IR divergences only at higher loop orders than considered 
here (this is checked explicitly in sec. [8]). Therefore it is guaranteed, from general principles, 
that all IR divergences must cancel in the spectral function, and we are allowed to handle 
them through dimensional regularization, as long as it is applied consistently. Of course, 
the cancellation of IR divergences can also be checked directly [7J. At the current order, the 
cancellation in fact probably takes place in every "master" spectral function separately, and 
for an arbitrary value of /T, as has been observed in the fully bosonic case |28j . 

Having thereby concluded the discussion of Higgs corrections, we note that similar meth- 
ods work in all other well. The full list of structures appearing can be found in 
appendix A, and the corresponding spectral functions are given in appendix B. Whereas for 
the current paper the cuts of appendix B are sufficient, in some contexts it may be desirable 
to also know the ultraviolet expansions in the full Euclidean domain (cf. e.g. sec. [8]); these 
are listed in appendix C, and have been derived with methods explained in refs. |161 117j. 



''if the exponentially small terms are included, then eqs. (|5.6|l . (|5.9|) . (|5.14|l imply (for k = 0) 

^i-bPAbT^ r, , ,ko ko ko r , .ko , Jcq, 



Im 



1 = ^ |l + -B(y)-np(-)--[nB(y) + np(-)J|. (5.16) 
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6. Next-to-leading order analysis 



Including now also contributions from the top quark and from gauge bosons, we turn to 
our full NLO expressions. In order to obtain a "universal" representation, we make use of 
completions of squares and substitutions of integration variables in order to express the results 
in terms of a minimal number of independent "master" sum- integrals, listed in appendix A. 
We specify graph-by-graph results in NDR in t«,ns of these r^lstersg Gauge parameter 
independence (with respect to both gauge groups) has been checked separately, so here only 
the Feynman gauge results are shown. They read 



• ) 



d ' 



(6.1) 



O 



2|/i;/Bp|^tBp^c( 2Xb — 2Xc — 2Xd +2^c — + 



^!^Bp(5iB +3c/2b) 



(6.2) 
(6.3) 



I^I^Bp(fi'lB +^92b) 



Xu — X„ — X, 



X,+Xf 



(6.4) 



rD-2 



(6.5) 



K^WoIb + ^g'L) i^e - + Xg - Xh - 2Xh + Xj j , 



(6.6) 



with dashed, solid, doubled, and wiggly lines representing scalars, leptons, quarks, and gauge 
bosons, respectively. 

Inserting the cuts, or spectral functions, from appendix B; setting D = 4 — 2e; and renor- 
malizing according to 



\h r 

\ii-ub\ 



hy{ii)\ ii" Zy , with 



1 + 



(47r)2e 



\ht?N,--{g\^'igl)\ + 0(/) 



(6.7) 
(6.8) 



*We have checked that the same results are obtained, for every diagram, with the recipe described below 
eq. (|3.6p . whereas in the strict 't Hooft - Veltman scheme there are additional terms; cf. appendix D. 
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where denotes a generic renormalized coupling, we obtain 



Svr 



ImnE(/C) 



12A r n 



1 



(4^) 



1 + 



/C6 



+ (5? + 3gi 



/i2 29 \ kl + /cV3 /■ p (17np - IQn^) 
4(47r)2 V ^ ^ y j ^ ^ 



(6.9) 



We note in passing that the structure + fc^/S originates from terms hke the last one in 
eq. (|5.5p . which after averaging over the directions of p yields 



.3-2e 

A physical interpretation for this structure is given after eq. (|8.8p . 



(6.10) 



7. Next-to-leading order results 



If we choose the renormalization scale as /i = M in eq. (j6.9p ; denote the corresponding 
renormalized coupling by set A'c = 3; and insert eq. ()5.2p together with 



2t^4 



"so" 



240 



(7.1) 



then the "decay rate" from eq. (|2.3p becomes 



r(/c) 



1 



87r^/FTM2 

+ (5? + 3(7i) 



Ar2 

M2 
29 
8(47r) 



' 80 



21 



2(47r) 



+ 



7^2 / r4 4 A;2r4 



60 



4 

M4 + 3 M6 



M4 ^ 3 



M6 



+ O 5' 



4 



9 



M2 ' M6 



(7.2) 



This constitutes our main result. (The 0{g^) correction is worked out in sec. [8] and amounts 
to Ar2 \T{T - Smn/vr), where mn is the thermal mass parameter given in eq. (j8.5p .) 

Equation (j7.2p can be compared with ref. [7]. We completely agree on all T- independent 
terms, as well as on the fact that gauge corrections lead to no term proportional to T"^. As 
far as the Higgs contribution is concerned, we find a result larger by a factor 2. As far as the 
top correction is concerned, we find that all terms of 0{T'^) cancel like in gauge corrections, 
whereas in ref. [7j only a partial cancellation was observed, so that a term of 0{T'^) was left 
over. The corrections of 0{T'^) in eq. (17. 2p were not considered in ref. [7J. 
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As shown in sec. O we have cross-checked the Higgs contribution by independently comput- 
ing the fuh rate for a general vrT/M and taking the non-relativistic limit only afterwards, cf. 
eq. (j5.16p . As far as we can judge, the Higgs correction in ref. [7] was inferred from the effect 
that a thermal mass would have on a vacuum- like result; but since there are no IR issues at 
NLO, results emerge from momenta p ^> gT, and must be insensitive to thermal masses. In 
fact, taking a zero-temperature decay rate ~ — mg, with mg denoting a vacuum mass, 
and denoting by m| = the thermal mass, then sees. 15.21 and 15.31 show that there are two 
contributions, amounting to — uiq — )• M"^ — mg — m| — m|^j^mg, rather than a single 

as was assumed in ref. [7j . One of the contributions comes 
from a "double-pole" cut and requires a careful analysis. 

Let us finally consider the total rather than the differential production rate, obtained 
according to eq. (12. 3|) as 



mass shift exp(m|-p-3-)(M 



dm?. 



1{T) 



dN_ 



d^k 

(2^)3 



2nF(A:°) r(/C) 



(7.3) 



A straightforward integration with T from eq. (j7.2p yields 

7(T) : 



\hJ'^M^T(, . /M 



3c2T /M 



o 



, rp ,1 2M 



(7.4) 



with the coefficients 



\ht\ 
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2(47r) 



+ 



7^2r4 



60M4 



+ {gl + 35i) 



60M4 



29 



8(47r) 



80M4 



(7.5) 
(7.6) 



Uncertainties are like in eq. (17. 2p . except that C2 has only thermal corrections. Embedding 
eq. (]7.4p in cosmology, the number density, n(T) = ^3^, conventionally normalized to the 
total entropy density, s(T), evolves aj^ 



^ d /n(T) 
dT V s(r) 



7(T) 



3c2(r)s(T)i7(r) 



(7.7) 



where H{T) is the Hubble parameter. This equation is correct as long as the differential 
abundance remains below its equilibrium value at each k separately, cf. e.g. ref. [19]; otherwise 
a "back reaction" from Pauli exclusion needs to be included. 



^ The speed of sound squared, c^, is often approximated as i , which is true in a conformally symmetric 
situation, but not when mass thresholds are crossed or effects from running couphngs are taken into account. 
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8. OPE representation 



In sec. [5l it was mentioned that if the EucUdean correlator 11^ is considered, then some of 
the individual terms appear to be IR divergent. As pointed out in ref. [15], the nature of 
these divergences can be understood by representing the result in an OPE form. 

We note first that, carrying out a naive dimensionally regularized NLO computation, the 
(bare) condensate {4>^4>)rp reads 



2 - 12 Ab X, + 2 1 /i^B I ^iVc (2 Xc -le) + (5?B + SffL : 



1 - D 



Here appear the same IR divergent structures as were encountered in sec. [5j If we now write 



= + n^'', where 11^'' denotes the vacuum part, then by inserting the expansions 
from appendix C into eqs. (16.1 
part can be expressed as 



AT) 



t(0) 



6.6p and renormalizing according to eq. ()6.7p . the thermal 



1 



2e 



naive 
T 



where 



1 



3A 

3^ 



1 + 



(47r)2e 



6A + \ht\ lyc 



Nr--{gl + 3g^,) 



+ (8.3) 



happens to be the renormalization factor related to the Higgs mass parameter, m| = rriQZm- 
The prefactor in eq. ()8.2p is IR-safe (temperature independent) as is typical of a Wilson 
coefficient; all IR-sensitive terms of 0(T'^) have been eaten up by the condensate. 

In order to compute the condensate correctly at 0{g^) also in the IR, we need to add 
to eq. (|8.ip . which in this naive form represents an ultraviolet contribution from "hard" 
momenta p > vrT, the contribution from "soft" momenta, in particular from Matsubara zero 
modes. This type of computations are best formulated within dimensionally reduced effective 
field theories. In fact the result can be extracted e.g. from ref. [29], and reads 



27r 



(4^) 



6A + 



E 



1 



+ 



+ 0{g' 



(8.4) 



where the Standard Model thermal masses read, for Nc = 3, 



TTT-t; 



ml + 



A , \ht\ 



+ 



16 



11 



rrir 



glT^ 



11 



rrir 



2rp2 



-g2T 



(8.5) 
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Summing together eqs. (j8.1|) and (j8.4p . multiplying by Zm as required by eq. (|8.2p . and for 
convenience also resumming the Debye mass contributions of 0{g^) into m|, we get 



rp2 rp2 
+ 



E 



+ 



487r2 
3(5? + 3^1 



6A 



In 



In 



levrT 

7e \ 

47rr y 

4^/~3 



|/.,PiV.lnf^'l 



27e 



8ttT J 
2C'(-1) 



C(-i) 



+ 41n 



'27rr 



(8.6) 



Apart from eqs. (|5.2p . (j5.3p . the fermionic /p 



nF(p) _ T^M" 



■{l + 26[ln3^ + l + g^]} 



p 24 I" ' -^L""" SttT ' " ' C(-l) . 

needed here. All ultraviolet divergences have nicely cancelled out in eq. ()8.6p . 



was 



The thermal correction to the "decay rate" of eq. (|2.3p now comes from the cut of the 
Wilson coefficient in eq. (j8.2p . It is given by 

[A + Oig')] Zm{ct>^<p)T + 0{\K\^g^T^) . 



r(^)(x:) 



3|/^^ 



47r\/FTM2 



(8.7) 



Looking back at eq. (j8.6p . we observe that there is a subleading 0{g^) contribution to the 
decay rate (from the term — "^'^'^ ~ C>{g) in Zjn{4>^ (t>) t) which is sensitive to thermal resum- 
mations, however in a computable way. 

Let us finally briefly remark on the terms of O(T^) in eq. (18. 7p . As was pointed out in 



ref. [15], these contain expectation values of various components of the energy-momentum 
tensor. At finite temperatures, the most important contribution turns out to emerge from 
its traceless part: denoting {Q^'^)x = diag(e,p, p,p), we may define 

1 



fJ.U\ 



3 



1 



a/T 



diag{e, p,p,p) - - diag(H )(e - 3p) 



^^e + p) diag(l, i , i , i ) . 



(8.8) 



Since e + p = Ts, where s is the entropy density, this condensate vanishes at zero temper- 
ature. The non-trivial momentum dependence in eq. (j6.9p comes from terms of the type 
Y^- K.^K,p{Q^^)j,/lC^, just like in ref. [15], where denote various gauge-invariant subparts 
of the full tensor (we have not worked out the decomposition because it is not needed here). 
Thermal resummations affect e and p first at 0{g^T^) [30] , and thus r(^)(/C) at 0{g^T^). 



9. Conclusions 

The purpose of this paper has been to compute the thermal production rate of right-handed 
neutrinos in the non-relativistic regime, meaning at temperatures much below their mass. 
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but still higher than the electroweak scale: mtop ^ vrT <^ M. (Equivalently, one can speak 
of a decay rate, cf. eq. ()2.4p .) In this regime, the results can be organized in the form of 
an Operator Product Expansion, with successively higher powers of vrT representing thermal 
expectation values of gauge-invariant condensates [15]. We have worked at NLO in the 
coupling constants (even at NNLO in sec. [8]), and to the third order in an expansion in 
(ttT/M)^; the main result is shown in eq. (17. 2p . At order (ttT/M)^ we fully confirm previous 
NLO results in the literature [7j, if employing naive dimensional regularization for handling 75; 
at order (ttT/M)^ we confirm the absence of corrections from gauge bosons, but find different 
results for corrections from the Higgs scalar and the top quark. Our results of (vrT/M)^ are 
new, and display the leading non-trivial dependence on the spatial momentum (the terms 
proportional to in eq. (\7.2h . which contribute at {ttT /M)^ to the total production rate). 

Despite small differences, the numerical magnitude of our corrections is similar to what 
was found in ref. [7]. Even though our NLO thermal Higgs correction is larger by a factor 2, 
we also find an NNLO term, cf. eq. ()8.6p . which numerically cancels about half of the Higgs 
correction. The most substantial difference is that we find a complete cancellation of top 
corrections at ©(ttT/M)^. Therefore our thermal corrections are negative, cf. eqs. (17. 2p and 
()8.7p . not positive as was found for the top correction in ref. [7]. Their precise numerical 
infiuence on leptogenesis computations remains to be inspected but probably the effects are 
no larger than was found in ref. [7] because the largest (top) term is absent at ©(vrT/M)^. 

In scenarios of TeV scale leptogenesis, the rate that we have computed could be relevant 
down to temperatures around the electroweak scale. Remarkably, eq. (|8.7p shows that then 
the rate could be extracted non-perturbatively from Euclidean lattice simulations of the type 
that were developed for studying the electroweak phase transition in the 1990s. In addition 
eq. ()8.7p shows that terms of O^irT/M)'^ are necessarily proportional to A; this can be traced 
back to a mismatch of the renormalization factors of the neutrino Yukawa coupling (eq. (|6.8p ) 
and the Higgs mass parameter (eq. ()8.3p ) which ultimately leads to the logarithm in eq. ()8.2p . 

As an outlook, we envisage that taking the representations in eqs. (I6.ip - (l6.6p as starting 
points, it is a feasible if hard task to extend the NLO analysis to the relativistic regime, 
vrT ~ M. Indeed a similar step has previously been taken in QCD, by going from OPE-type 
results in ref. [K] to relativistic results in ref. |28] . We believe that in the current context it 
would be important to keep the spatial momentum k different from zero, which makes the 
analysis more demanding; still, in principle similar techniques should work. 
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Note added 



In the revised version (v3) and Erratum of ref. [7j, the results have been corrected and now 
agree with ours. 



Appendix A. Definitions of master sum-integrals 

Denoting by P, Q bosonic and by K, R, S fermionic Matsubara four-momenta, and employing 
the usual conventions and ^ for the corresponding measures, the master sum-integrals 
entering the computation are defined as follows (in the graphical notation, a dashed line 
indicates a bosonic propagator, a solid line a fermionic one, a filled blob a squared propagator, 
and a cross that the momentum appears in the numerator as well): 



J. 



J. 



1 

P 

1 

m 



PQ 



Q^P^ 
PQ 

1 



R'^P^ 
P{R} 

1 



P^R-i 
P{R} 

1 



' f^^^^R^P\P-K) 



A.l 
A.2 

A.3 
A.4 
A.5 
A.6 
A.7 
A.8 
A.9 



{RS} 



16 



^2 



P{R} ^ ' 



I-----H Xe 



{RS} ^ ' 



^^^^R^{R- Pf{P - Kf ' 
^ t^^Q^P^{Q-PY{P-KY ' 



K 

P{R} 



X. 



P{R} 



P{R} 

X 



R'^P'^{R 


-PY{P- 










R'^P^iR 


-py{R- 


KY ' 




2K-P 




R'^P^iR 


-py{R- 











Appendix B. Results for master spectral functions 

The spectral functions corresponding to the structures of appendix A are obtained from 

PXy, = Im[X^];.,^^_i[;i;0+io+] • (B.l) 

As has been elaborated upon in sec. [5l they can be derived by carrying out the Matsubara 
sums, which corresponds to cutting zero, one, or two lines and weighting them with a thermal 
distribution. In the result, cuts correspond to logarithmic terms, which can only arise from 
vacuum loops. Given that 2-cut contributions have no vacuum loop left, spectral functions 
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only arise from 0-cut and 1-cut contributions. After an expansion in small thermal momenta, 
the vacuum parts always have simple cuts, and thereby the spectral functions corresponding 
to the master structures of appendix A are elementary, cf. eqs. ()B.4p - ()B.2ip below. 

Before presenting the list, we would like to mention that apart from the procedure described 
in sec. O there exists also a more elaborate possibility for determining the spectral functions. 
Namely, one can first "blindly" work out the OPE expansion in the Euclidean domain up to 
the desired order in l/K"^ . The corresponding results are listed in appendix C. In these results, 
non-analytic K-dependence only appears in a factor X, defined in eq. (IC.2p . Introducing the 
MS scale parameter, it can be expanded as 



X 



1^ 



'2e 



(47r)2 



1 



12 



and the corresponding spectral function then arises from 

,-2e 



Im[^]A;„^-i[feo+iO+] = sign(A;°) ^ 



167r 



(B.2) 



(B.3) 



So, we observe that spectral functions can only arise from structures of type X"/e™, and 
that all Euclidean structures that are either finite or have 1/e-poles with the non-analytic 
scale dependence given by T rather than by K, yield vanishing spectral functions. 

Proceeding to the list, we wish to remove clutter by not showing the non-consequential 
factor /U"^*^ explicitly, and by similarly omitting the arguments of the functions n^, Up. In 
addition, the errors of the 2-loop structures, which are 0{e, ^), are not displayed. Thereby 
the spectral functions corresponding to the master structures read 



PXc 
P^A 
Pt. 



pj. = 0' 



167r 
p IGvrp 



P 



l + e( ln^ + 2 



p IGvrp ' 



P% = Pic 



rip 



p 16vrp ' 

Up 

p 167rp ' 



p IGvrp ' 

rip 



(B.4) 
(B.5) 
(B.6) 

(B.7) 
(B.8) 
(B.9) 

(B.IO) 

(B.ll) 

(B.12) 
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= p~^= (B.13) 
K? f 2nB -np , . 

_ K? f 3np 

/■ np-2nB , A;g + A;V3 /" /Unp TneN 



Appendix C. Euclidean large- momentum expansions 

For completeness, we list here Euclidean large-momentum expansions for the master sum- 
integrals defined in appendix A. In order to be as concise as possible, we introduce the 
notation 

{n-B , bosonic Hne -, n 

— np , fermionic line 

and provide expressions valid simultaneously for all the statistics carried by the lines. The 
shorthands 

K-^^ r(i + 6)r2(i-6) r(i + 26)r2(i-26) 

^ (47r)2-^ r(l-2e) ' ~ r(l - 3e)r2(l + e)r(l - e) ' ^ ' 

as well as 



Qk = " (C.3) 
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are also helpful. Thereby we obtain 



e ^ e 



p 



O 



e(l-2e)jp p 
l-2e /■ ni{p)n2iq) 



1 ^ 



1 -^p.q 



p,q 



+ 



+ 



n3{p)ni{q) 

'p P ^ Jp,q 

e + 2(l-e)(3-2e)ei^ f n3(p)ni((7) 



(2 - e)i^2 
1 /■ n3{p)n2{q) 



p,q 



p,q P9 



2e(l -2e)(l -3e)(2-3e) 



+ 



+ 



— t 



ni(p) 2e(l + e) GA'pni(p) 



+ 



E 



ni{p)nj{q) 
pq 



+ 



ji6 



i^2x2 

£2(1 -2e)^ 



+ 



+ 



X 



E 



p 



_}_ f {ni+n2){p){ns+n^){q) / 
K^L^ pq 
K^X^Y 



2e2(l -2e)(l -3e) 



+ 



X 



e(l - 2e) 7p 



ni(p) , 2il + e){2 + e)eKpni{py 



P 



+ 



i^2 



E 



=3 •'P 



ni(p) , 2{l + e)i2 + €)QKPni{p) 



+ 



3i^2 
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+ 



+ 



+ 



+ 



+ 



+ 



+ 



+ 



e(2 - e)K^ J^^^ pq 



ni{p)ni{q) .^/^S 



2e2(l -3e)(2-3e) 
X 



ni{p) ^ 2{l + efQKpni{jp) 



P 



X f 2{l + e)eKpn3{p) 



2X 
1 -2e 



n^jp) ^ (1 + 6) (1 + 2e2) pn4(p) 



p 3eK^ 
e + 2(1 - e)^ 9^- f (n2 + n3)(p)n4(g) - n2{p)n-i{q) 

6(2 - 6)^:2 y^^^ 



pq 



2 



ni(p)^4(g) 



2 v2/ 



j^^x"(i-y) 

e3(l - 2e) 
2X(l + e) 



ni{p) ^ 2(ll + 6e + e2)ei^pni(p) 



3K2 



n5(p) ^ (2 + e )(3 + e) Qxpn^ip) 



p 3i^2 
e + 2(1 - e)2 9j^ /" ni{p)n2{q) + n2.{p)ni{q) 

6(2 - 



e(l + e) + 6(l-e)2ei^ 



e(2 - e)E:2 



i=i -^P'q 



P9 

ni{p)n^{q) 
pq 



1 



pq 



^i(p)"3(g) +n2{p)ni{q) ^ c(^^ ^ 



pq 



(C.12) 



(C.13) 



(C.14) 



Appendix D. On the treatment of Dirac traces 



In analogy with the leading-order example of eq. (j4.6p . let us compare the NDR expressions 
in eqs. (|6.ip - (|6.6p with the recipe explained below eq. (j3.6p as well as with the strict 't Hooft - 
Veltman scheme. After Wick contractions, Lorentz algebra, isospin traces and using standard 
properties of non-chiral Dirac matrices (but doing nothing with 75 for the moment), the 
contributions to the Euclidean correlator IiE{K) in Feynman gauge read (omitting the overall 
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factor l/ijyBp) and abbreviating Lorentz-indice^ /ii, ...,/X6 as in 71 = 7^^ or K2 = K^^ etc.) 

12ABTr[7ia,72aJ^ , (D.l) 

2|/iiBpiVc Tr [7ia,72 a,] TV [73a,74 a,] ^ ^^^5 " ol'S.^o " m2 ' (D-2) 




„ ^^^^ P4(^ _ p)2R2^p _ py 

2 , o 2 ^^ r K^{K-P)2 



= (SiB + 352b) Tr [710172 Or] y y Q2p4(^K - Py ' ^^'^'^ 

,<->. /2^^2^^.^ . 1 vf i^l (i^ - P)2 (P + Q)^ 

= -(5lB + 352b) Tr [710^72 Or] - p4Q2(p_Q)2(j^_p)2 ' (^-4) 

^ Iv/; i^l(j^-P)2(i^-Q)3(i^-P)4 

S1234 2 jL^^p2(p_Q)2(^_p)4(^_Q)2 ' ^^-^^ 

S1234 = (5iB + 352b) - 2) Tr [710^76 (523546 - 524536 + 534526) 

- 5iB Tr [710^7276 aR7376 aR74aR] 

- 352BTr [710^7276 (aL7376aL- 7376) 7401,], (D.6) 
, 1^ Kr{K-QUP + QUK-P), 

11234 2|.^^p2Q2(p_g)2(^_p)2(^_g)2 ' ^ ''l 

ii234 = (5iB + 352b) Tr [71 aL76 Or] (523546 - 524536 + 534526) 

- 352b Tr [7iaL7273aR74aR] ■ (D.8) 

All diagrams with Higgs self-energy insertions are seen to be proportional to the same struc- 
ture that already appeared in the LO contribution; if handled as shown below eq. (14. 6p . they 
lead to Tr [7iaL72aR] 2512 and immediately reduce to the NDR results. On the other hand, 
the two diagram classes with the vectors coupling to the lepton line do seem to get additional 
contributions proportional to different Dirac traces, which have been separated in the second 
and third lines of eqs. (1D.6P and (1D.8|) . 

It turns out, however, that upon employing the prescription explained below eq. (|3.6p these 
additional contributions vanish identically, which leaves us with the same results as in NDR 
for each of the diagrams. As a specific example, consider the second line of eq. (jD.Sp . 

8 Tr [7i«L7273aR74aR] = Tr [717273(74 + 757475)] - Tr [71757273(74 + 757475)] 

- Tr [71757273(7574 + 7475)] + Tr [717273(7574 + 7475)] 
= Tr [7471727375 + 7172737475] , (D.9) 

where we have used an anticommuting 75 as well as 7I = 1 in traces with more than one 75. 
Then, for the remaining traces with a single 75 we get 4(64123 + ^1234), such that the two 



^Note that we do not use /is here, to avoid confusion of 7;j5 with 75. 
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terms cancel due to antisymmetry of the e tensor. In a completely analogous way, the last 
two lines of eq. ()D.6P are seen not to contribute in this specific scheme. 

If the 't Hooft - Veltman scheme is used rather than the recipe below eq. ()3.6p . then graphs 
with closed fermion loops do differ from those in NDR. To see this, note that the integral in 
eq. (|D.2p evaluates at zero temperature to 



KiiK 
PR P\K 



P)^R^{P-Ry 



144 



1 31 763 , 



4e(l - 2e)^^^|^ + (5 + 2e)K2 534 - (1 - 2e)(i^3 <724 + K^923) 



(D.IO) 



where the coefficients X and Y are defined in eq. (|C.2p . With Tr [730^74 ^r] — ^ 2^34 this 
yields the NDR result, but with eq. (|3.3p we rather have Tr [73aL740R] = 2534, where g^i, = 1, 
if yu = < 3, and Qf^u = otherwise. Contracting with the 334 in eq. (ID.lOh produces 4 
rather than 4 — 2e, which turns into a difference of 0{1) in the spectral function because of 
the prefactor 1/e^ (concretely, | in eq. (16. 9p turns into y|). That said, this difference can 
presumably be "hidden" if the neutrino Yukawa coupling is expressed in terms of a physical 
quantity (such as a pole mass) through a computation carried out in the same scheme. 

In thermal corrections to eq. (jP.lOp . there is only a prefactor 1/e in the Euclidean do- 
main (cf. appendix C), which implies that the corresponding spectral function is finite (cf. 
appendix B). Therefore we expect that the ambiguity of 0(e) in a prefactor does not affect 
NLO thermal corrections to spectral functions. Nevertheless, it might be interesting to work 
out the full tensor integrals of eqs. (lD.ip -( lD.7p in the OPE regime, in analogy with the ex- 
pansions in appendix C, thereby producing results for any desirable scheme. Unfortunately 
this involves a substantial amount of work and goes beyond the scope of the present study. 
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